%% 统计学预测方法Statistical data prediction methods
%% 1.典型微分方程、指数平滑法和灰色预测
% Typical differential equation, exponential smoothing method and grey prediction
%%  1.1.典型微分方程
% 微分方程是数学建模的重要方法，许多实际问题的数学描述将导致微分方程的定解问
% 题。列方程常见的类型有3种：(1)按规律直接列方程；(2)微元分析与区域积分；(3)
% 模拟近似法。
% 在各类学科都有很多经过实践检验的规律和定律可以列写微分方程，如牛顿运动定律、
% 基尔霍夫电流电压定律、物质放射性规律、曲线切线性质等。
% 例(物体冷却过程)：牛顿冷却定律指出，当物体表面与周围存在温度差时，单位时间
% 内从单位面积散失的热量与温度差成正比，比例系数为热传递系数，记为k，试建立
% 相关传热模型。
% (2018数学建模国赛A.1)高温作业服可避免人们在高温环境下作业受到灼伤。防护服分
% 为4层，其中Ⅰ、Ⅱ、Ⅲ层为织物制造，第Ⅳ层为织物与皮肤之间的空气间隙；第Ⅰ
% 层直接与外部环境接触,假人体内恒温为 37℃，外界温度为75℃。试建立非稳态传热
% 模型，反应不同时间节点的传热情况。
% 数据见文件:data1。
% 利用文件中附件1信息可以建立逻辑严密的热传导模型，读者不妨思考，可不可以只
% 用附件1有关厚度信息，忽略边界效应直接利用牛顿冷却建立微分方程，利用优化算法
% 去逼近附件2中数据呢？答案是可行的。
% 设最外层防护服与空气间热传递系数为k0，Ⅰ、Ⅱ、Ⅲ、Ⅳ层内部热传递系数分别为
% k1,k2,k3,k4，Ⅳ与皮肤间热传递系数为k5，问题变为建立模型，求Ⅳ靠近皮肤的表面
% 温度，使其值逼近附件2。
% 可编写如下程序建立简单的传热模型：
% (该程序可见文件：Newton_cooling_case.m)
clear,clc
k0=1;k1=0.99;k2=0.98;k3=0.9;k4=0.99;k5=0.01;% 传热系数
d=[0.6,0.6,3.6,0.6];% 厚度设置
d1=cumsum(d);
dx=[1e-2,1e-2,1e-1,1e-2];% 厚度微元设置
tw=75;% 外部温度75
tn=37;% 假人温度37
u1=tn*ones(1,length(0:dx(1):d1(1)));
u2=tn*ones(1,length(d1(1):dx(2):d1(2)));
u3=tn*ones(1,length(d1(2):dx(3):d1(3)));
u4=tn*ones(1,length(d1(3):dx(4):d1(4)));% 防护服温度初始化
tm=5400;% 总时间
u=zeros(1,tm);
for t=1:tm
    tic
    for i=1:length(0:dx(1):d1(1)) % 第一层防护服温度迭代
        if(i==1)
            deltau=-k0*(u1(i)-tw)-k1*(u1(i)-u1(i+1));
            u1(i)=u1(i)+deltau;
        elseif(i==length(0:dx(1):d1(1)))
            deltau=-k1*(u1(i)-u1(i-1))-k2*(u1(i)-u2(1));
            u1(i)=u1(i)+deltau;
        else
            deltau=-k1*(u1(i)-u1(i-1))-k1*(u1(i)-u1(i+1));
            u1(i)=u1(i)+deltau;
        end
    end
    for i=1:length(d1(1):dx(2):d1(2)) % 第二层防护服温度迭代
        if(i==1)
            deltau=-k1*(u2(i)-u1(end))-k2*(u2(i)-u2(i+1));
            u2(i)=u2(i)+deltau;
        elseif(i==length(d1(1):dx(2):d1(2)))
            deltau=-k2*(u2(i)-u2(i-1))-k3*(u2(i)-u3(1));
            u2(i)=u2(i)+deltau;
        else
            deltau=-k2*(u2(i)-u2(i-1))-k2*(u2(i)-u2(i+1));
            u2(i)=u2(i)+deltau;
        end
    end
    for i=1:length(d1(2):dx(3):d1(3)) % 第三层防护服温度迭代
        if(i==1)
            deltau=-k2*(u3(i)-u2(end))-k3*(u3(i)-u3(i+1));
            u3(i)=u3(i)+deltau;
        elseif(i==length(d1(2):dx(3):d1(3)))
            deltau=-k3*(u3(i)-u3(i-1))-k4*(u3(i)-u4(1));
            u3(i)=u3(i)+deltau;
        else
            deltau=-k3*(u3(i)-u3(i-1))-k3*(u3(i)-u3(i+1));
            u3(i)=u3(i)+deltau;
        end
    end
    for i=1:length(d1(3):dx(4):d1(4)) % 第四层防护服温度迭代
        if(i==1)
            deltau=-k3*(u4(i)-u3(end))-k4*(u4(i)-u4(i+1));
            u4(i)=u4(i)+deltau;
        elseif(i==length(d1(3):dx(4):d1(4)))
            deltau=-k4*(u4(i)-u4(i-1))-k5*(u4(i)-tn);
            u4(i)=u4(i)+deltau;
        else
            deltau=-k4*(u4(i)-u4(i-1))-k4*(u4(i)-u4(i+1));
            u4(i)=u4(i)+deltau;
        end
    end
    toc
    u(t)=u4(end);% 第四层靠近皮肤表面温度
end
% 此时需要做的工作是利用最小二乘编写优化算法，优化各层传热系数，使其逼近附件
% 2值即可。
% 例(导弹追击问题)：设为一导弹发射点A(0,0),(0,10)处有目标物向右以恒定速率v运
% 动，0时刻处发射一枚导弹以2v的速度追赶目标物，试通过编程，输入b、v，输出导
% 弹动态路线。
clear,clc
b=input('请输入初始时刻目标物的高度:');
v=input('请输入目标物运动速度大小:');
m=0;% 目标物横坐标
q=0;% 导弹横坐标
p=0;% 导弹纵坐标
dt=1;% 时间微元
for i=1:1000
    plot(m,b,'b*-',p,q,'rx--');% 绘制当前导弹和目标物位置
    axis([0,p+10,0,b+1]);
    hold on
    pause(0.2);% 动态显示
    m=m+v*dt;% 目标物横坐标运动
    d=sqrt((m-p)^2+(b-q)^2);% 目标物与导弹间距
    p=p+(2*v/d)*(m-p)*dt;% 导弹横坐标运动
    q=q+(2*v/d)*(b-q)*dt;% 导弹纵坐标运动
    if d<0.05% 当两者距离小于某值时停止
        break
    end
end
hold off
% 例(传染病SIR模型)：假设人群分为健康者，病人，和病愈后具有免疫力不得病三类
% 设任意时刻t，这三类人群占总人口比例分别为s(t),i(t),r(t)，设病人的日接触率
% 为lamda，日治愈率为miu，则有传染强度delta=lamda/miu，若人口总数为N，则有微
% 分方程组：
% s(t)+i(t)+r(t)=1;dr/dt=miu*i;ds/dt=-lamda*s*i;di/dt=lamda*s*i-miu*i;
% 利用matlab仿真此类传染病的发病过程。
clear,clc
N=1e6;% 人口总数
s0=0.98;% 健康者初始比例
i0=0.02;% 病人初始比例
r0=0;% 免疫者初始比例
dt=1;% 单位时间/天
lamda=3;% 病人日接触率
miu=0.5;% 日治愈率
s=s0;i=i0;r=r0;
for k=1:1e2 % 总天数
    plot(k,N*s,'b*-',k,N*i,'rx--',k,N*r,'bx--');% 绘制当前各类人群人数
    axis([0,k+10,0,N]);
    legend('健康者','病人','免疫者');
    hold on
    pause(0.2);% 动态显示
    i=i+(lamda*s*i-miu*i)*dt;
    s=s+(-lamda*s*i)*dt;
    r=r+(miu*i)*dt;
end
hold off
% 利用该模型可以简单模拟传染病发病过程。
% 对于常微分方程，matlab提供了几个解常微分方程数值解的函数，如：ode45，ode23
% ，ode113等，其中ode45采用四五阶龙格-库塔法(RK方法)，是解非刚性微分方程的首
% 选方法，ode23采用23阶RK方法，ode113采用多步法，效率一般比ode45高。
% 例(Logistic-弱肉强食模型)：对于一种群，当存在自然资源，环境条件限制时，种
% 群增长率将会随种群数量增加而减小，设t时刻种群数量为x(t)，种群增长率r(x)为x
% 的线性函数，即：r(x)=R-s*x;则不考虑其他条件时，该种群数量增长可表示为：
% dx/dt=r(x)*x;
% 当自然界存在捕食者和被捕食者两种群时，设t时刻两种群数量分别为x1(t),x2(t)，
% 种群增长率分别为r1(x1)=(1e2-x1)/(50),r2(x2)=(1e3-x2)/(100);,捕杀量为
% 0.2*x1*x2;利用matlab仿真两种群数量变化。
clear,clc
x10=5;% 捕食者初值
x20=500;% 被捕食者初值
dx=@(t,x)[(1e2-x(1))/(50)*x(1)+0.01*x(1)*x(2);
    (1e3-x(2))/(100)*x(2)-0.005*x(1)*x(2)];% 列写微分方程
[t,x]=ode45(dx,[0,100],[x10,x20]);
plot(t,x(:,1),'b*-',t,x(:,2),'rx--')
legend('捕食者','被捕食者')
%%  1.2.指数平滑法
% 本节及之后讨论对时间序列数据的预测问题。假设预测时，离预测点越近的点作用越
% 大，随着时间变化权重以指数形式变化，则可以采用指数平滑法。一般情况下，时间
% 序列可以分解为水平部分(均值)，趋势部分(上升或下降)，季节性部分(周期重复)，
% 随机波动(白噪声)。
% 当数据没有趋势和季节性时，可以采用一次指数平滑法；当数据有趋势无季节性时，
% 可以采用二次指数平滑法；当数据既有趋势又有季节性时，采用三次指数平滑法。
% 对于数据噪声，可以采用信号处理相关方法，如小波变换进行滤除。
% 现做以下规定：y(t)为t时刻数据，y'(t)为t时刻预测值，sn(t)为t时刻n次指数平滑
% 值，alpha为记忆衰减因子。
% 一次指数平滑：y'(t+1)=alpha*y(t)+(1-alpha)*y'(t);s1(t)=y'(t);
% 二次指数平滑：y'(t+T)=a(t)+b(t)*T;其中，a(t)=2*s1(t)-s2(t);
% b(t)=alpha/(1-alpha)*(s1(t)-s2(t));s2(t)=alpha*s1(t)+(1-alpha)*s2(t-1);
% 三次指数平滑：y'(t+T)=a(t)+b(t)*T+c(t)*T^2;其中，a(t)=3*s1(t)-3*s2(t)+s3(t);
% b(t)=alpha/(2*(1-alpha)^2)*((6-5*alpha)*s1(t)-2*(5-4*alpha)*s2(t)+(4-3*alpha)*s3(t));
% s3(t)=alpha*s2(t)+(1-alpha)*s3(t-1);
% alpha取值标准：水平趋势，0.05~0.2；有波动，长期趋势变化不大，0.1~0.4；波动
% 大，长期趋势变化大，0.6~0.8；上升或下降的发展趋势，0.6~1。
% 本节利用指数平滑法对具有上升趋势和周期性的数据进行预测：
% (程序可见文件：Cubic_exponential_smoothing.m)
clc,clear 
x=0:1e-1:20;y=2*exp(-0.16*x).*sin(2*x)+1/2*x;% 原始信号
alpha=0.6;% alpha值
n=length(y);% 原始信号长度
s10=mean(y(1:3));% 取y前三个数据均值作为各次指数平滑初值
s20=s10;
s30=s10; 
s1=alpha*y(1)+(1-alpha)*s10;% 各次指数平滑值初始化 
s2=alpha*s1(1)+(1-alpha)*s20; 
s3=alpha*s2(1)+(1-alpha)*s30; 
for i=2:n+20 % 预测原始信号后20位数据值
    if i>n
        y=[y,a+b+c];% T=1时，三次指数平滑预测值，此处进行滚动预测
    end
    s1=alpha*y(i)+(1-alpha)*s1;% 各次指数平滑值迭代
    s2=alpha*s1+(1-alpha)*s2;
    s3=alpha*s2+(1-alpha)*s3;
    a=3*s1-3*s2+s3;
    b=alpha/(2*(1-alpha)^2)*((6-5*alpha)*s1-2*(5-4*alpha)*s2+(4-3*alpha)*s3);
    c=alpha^2/(2*(1-alpha)^2)*(s1-2*s2+s3);
end
plot(1:n+20,y)% 绘制预测图
%%  1.3.灰色预测
% 灰色预测适合处理小样本，中短期，指数变化的数据。设有时间序列x(n),n=1,2,...,N
% ，计算序列级比lamda(k)=x(k-1)/x(k),k=2,3,...,N，若Any(对于任意) lamda(k) 
% Belongto(属于) (exp(-2/(N+1)),exp(2/(N+1)))，则可进行GM(1,1)预测。
% 现令:x1(n)=cumsum(x(n)),n=1,2,...,N为x(n)一次累加序列；
% x0(n)=diff(x(n)),n=1,2,...,N-1为x(n)一次累减序列：
% z(n)=1/2*(x1(2:end)+x1(1:end-1))为x1(n)均值序列。
% GM(1,1)预测：
% 求解：[a,b]'=(B'*B)\B'*x(1:N)';其中，B=[-z(1:N)',ones(N,1)];则可预测得到：
% x1'(k+1)=(x(1)-b/a)*exp(-a*k)+b/a;x'(k+1)=x1'(k+1)-x1'(k);
% 进行GM(1,1)相对误差检验：delta(k)=abs(x'(k)-x(k))/x(k),k=1,2,...,N,当
% delta(k)<0.1时，达到较高精度；delta(k)<0.2时，达到一般精度。进行GM(1,1)
% 级比偏差值检验：rho(k)=1-((1-1/2*a)/(1+1/2*a))*lamda(k);rho(k)<0.1时，达到
% 较高精度；rho(k)<0.2时，达到一般精度.
% GM(2,1)预测：
% 对于非单调摆动发展序列或有饱和s型序列，可以考虑建立GM(2,1),DGM,Verhulst模型：
% 称d2(x1(t))/xt2+a1*dx1(t)/dt+a2*x1(t)=b为GM(2,1)的白化方程。令：
% B=[-x(2:end)',-z(1:end)',ones(1:N-1)'];Y=[x0(1:end)'];
% 则有GM(2,1)模型最小二乘估计：[a1,a2,b]'=(B'*B)\B'*Y;
% 得到参数后利用边界条件x1(1),x1(end)求解白化方程即可进行预测。
% DGM(2,1)预测：
% 令B=[-x(2:end)',ones(1:N-1)'];Y=[x0(1:end)'];得到：[a,b]'=(B'*B)\B'*Y;
% 则可预测得到：x1'(k+1)=(b/a-x(1)/a)*exp(-a*k)+b/a*k+(1+a)/a*x(1)-b/a^2;
% Verhulst模型：
% 令B=[-z(1:end)',[z(1:end).^2]'];Y=[x(2:end)'];得到：[a,b]'=(B'*B)\B'*Y;
% x1'(k+1)=a*x(1)/(b*x(1)+(a-b*x(1))*exp(a*k));
% 本节给出matlab简要计算各灰色预测模型的代码：
% (代码可见文件:Grey_prediction.m)
% 灰色预测模型汇总
clear,clc
% 生成各类型序列
x=0:0.1:10;
f1=diff(exp(-0.02*x)+1);% % 指数下降类序列，适合GM(1,1)预测
f2=diff(1e2*exp(x)-1e1*exp(0.001*x)+2);% 双指数综合，可选用GM(2,1)预测
f3=diff(2.7*exp(-0.4*x)+4*x+1);% 指数与一次综合，可选用DGM(2,1)预测
f4=diff(1./(1+50*exp(-0.01*x)));% s型曲线，可选用Verhulst预测
f=[f1;f2;f3;f4];
% 预测数据，与原始数据作比较
for i=1:4
    data=f(i,:);
    n=length(data);% 原始数据长度
    data0=diff(data);
    data1=cumsum(data);
    z_data=1/2*(data1(2:end)+data1(1:end-1));
    if(i==1)% GM(1,1)
        B=[-z_data(1:end)',ones(n-1,1)];
        Y=[data(2:end)]';
        u=(B'*B)\B'*Y;
        pre=@(k)(data(1)-u(2)/u(1))*exp(-u(1)*(k-1))+u(2)/u(1);
        fp=pre(1:n);fp=diff(fp);
    end
    if(i==2)% GM(2,1)
        B=[-data(2:end)',-z_data(1:end)',ones(n-1,1)];
        Y=data0';
        u=(B'*B)\B'*Y;
        syms x_un(t)
        x_un=dsolve(diff(x_un,2)+u(1)*diff(x_un)+u(2)*x_un==u(3), ...
            x_un(0)==data1(1),x_un(n-1)==data1(end));% 求解白化方程
        xt=vpa(x_un,6);
        fp=double(subs(x_un,t,0:n-1));
    end
    if(i==3)% DGM(2,1)
        B=[-data(2:end)',ones(n-1,1)];
        Y=data0';
        u=(B'*B)\B'*Y;
        pre=@(k)(u(2)/u(1)^2-data(1)/u(1))*exp(-u(1)*(k-1))+...
        u(2)/u(1)*(k-1)+(1+u(1))/u(1)*data(1)-u(2)/u(1)^2;
        fp=pre(1:n);fp=diff(fp);
    end
    if(i==4)% Verhulst
        B=[-z_data',z_data'.^2];
        Y=data(2:end)';
        u=(B'*B)\B'*Y;
        pre=@(k)u(1)*data(1)./(u(2)*data(1)+(u(1)-u(2)*data(1))*exp(u(1)*(k-1)));
        fp=pre(1:n);fp=diff(fp);
    end
    figure(i)
    plot(data,'b*-')
    hold on
    plot(fp,'r*-')
    hold off
    legend('实际值','预测值')
end
% 从预测结果中可以清楚的观察到部分结果预测并不准确，灰色预测适合小样本数据，
% 当样本数据足够多时可选择其他时间序列预测方法，此外，样本变化幅度不能过大，
% 在预测前需要对数据进行检验，灰色预测只能进行中短期预测，长期预测结果将出现
% 较大误差。
%% 2.马尔可夫预测Markov predicted
% 某一系统在已知现在情况的条件下，系统未来时刻的情况只与现在的情况有关，而与
% 历史无直接关系，则该系统可称其位马尔可夫系统。
% 有马氏链模型：某序列x(n)满足：
% p{x(n+1)=j|x(n)=i(n),x(n-1)=i(n-1),...,x(1)=i(1)}=p{x(n+1)=j|x(n)=i(n)};
% 即未来时刻值只与现在值有关的序列。
% 时齐的马氏链：由一个状态转移到另一个状态的概率只与时间间隔有关，与起始时刻
% 无关。称以m步转移概率Pij(m)位元素的矩阵P(m)为马氏链的m步转移矩阵。
% 例：设某系统状态空间为：
E=[1,2,3,4];
% 有观察数据：
a=num2str([4,3,2,1,4,3,1,1,2,3,2,1,2,3,4,4,3,3,1,1,1,3,3,2,1,2,2,2,4,4,2,3,2,3,1,1,2,4,3,1]);
% 则有该马氏链一步转移矩阵P(1):
f=zeros(length(E));
for i=E
    for j=E
        f(i,j)=length(strfind(a,num2str([i j])));% 寻找a中出现[i,j]的次数
    end
end
p=sum(f,2);
for i=1:4
    f(i,:)=f(i,:)/p(i);
end% 一步转移矩阵
%% 3.ARIMA预测ARIMA Forecast
%%  3.1.平稳性Daniel检验
% 设有一组观察数据：(x(1),y(1)),(x(2),y(2)),...,(x(n),y(n))，设{x(1),x(2),...,x(n)}
% 的秩统计量为{R(1),R(2),...,R(n)}(秩统计量可理解为若将总数居从小到大排，该数
% 据排到第几),{y(1),y(2),...,y(n)}的秩统计量为{S(1),S(2),...,S(n)}。令：
% d={R(1)-S(1),R(2)-S(2),...,R(n)-S(n)}则观察数据有Spearman相关系数：
% Qxy=1-6/(n*(n^2-1))*sum(d.^2);
% 做假设检验:H0:rho(xy)=0,H1:rho(xy)~=0,其中rho(xy)为x,y总体相关系数。
% 令T=Qxy*sqrt(n-2)/sqrt(1-Qxy^2);，当abs(T)<=t(alpha/2)(n-2)时接受H0。
% 当H0检验通过时，说明观察数据序列平稳，不通过时说明存在上升或下降趋势。
clear,clc
% 设有数据列：
n=1000;
x=cumsum(ones(1,n));
y=randn(1,n);% 高斯噪声列
ry=sort(y);Ry=zeros(1,n);
for i=1:n
    [~,r]=find(y==ry(i));
    Ry(i)=r;% 获得y数据列秩统计量
end
d=x-Ry;
Qxy=1-6/(n*(n^2-1))*sum(d.^2);
T=Qxy*sqrt(n-2)/sqrt(1-Qxy^2);
t_0=tinv(0.975,n-2);% 计算上alpha/2分位数
% 由于abs(T)<t_0，则通过假设检验，数据列平稳。
%%  3.2.自相关系数与偏自相关系数
clear,clc
% 设有数据列：
n=1000;
x=cumsum(ones(1,n));
y=randn(1,n);% 高斯噪声列
y_mean=mean(y);% 数据y平均值
y_0=y-y_mean;% 原始序列减去其均值
% 则自相关系数可表示为：
rho=zeros(1,n);
for k=0:n-1
    rho(k+1)=sum(y_0(1:n-k).*y_0(k+1:n))/sum(y_0.^2);
end
% 当序列为噪声序列时，自相关系数rho(k+1)~N(1,1/n)。
% 偏自相关系数计算：
phi=zeros(1,n-1);
for k=0:n-2
    rho1=rho(2:k+2);
    kp=length(rho1);
    rho2=zeros(kp);
    for k0=1:kp
        rho2(k0,k0:end)=rho(1:kp-k0+1);
        if(k0>=2)
            rho2(k0,1:k0-1)=fliplr(rho1(1:k0-1));
        end
    end
    phi1=rho2\rho1';phi(k+1)=phi1(end);
end
%%  3.3.ARIMA
% AR(p)模型：p阶自回归，y(t)=miu+sum(r(i)*y(t-i))(i=1:p)+et,其中,et为残差，
% miu,r为待定系数。
% MA(q)模型：q阶滑动回归：y(t)=miu+sum(r(i)*e(t-i))(i=1:q)+et,其中,et为残差，
% miu为待定系数。
% ARIMA中I指差分，ARIMA模型处理的数据为平稳数据，当数据有上升下降趋势时，需
% 要通过差分获得平稳数据列，若令需要差分的次数为d，则有完整的ARIMA(p,q,d)模
% 型：y(t)=miu+sum(r(i)*y(t-i))(i=1:p)+sum(r(i)*e(t-i))(i=1:q)+et。
% 本节给出matlab ARIMA模型预测的简单方法：
% (程序可见文件：ARIMA_model_prediction.m)
% ARIMA模型预测ARIMA model prediction
clear,clc
% 原始数据列：
n=1000;
x=cumsum(ones(1,n));
y=1/n*x+sin(x/(4*pi))/10+randn(1,n);% 具有上升趋势，小幅周期波动,噪声的数据列
% 判断数据列平稳性并做差分处理
dm=4;% 约束最大进行4次差分
yd=y;% 差分数列初始化
for dp=1:dm+2
    n1=length(yd);
    ry=sort(yd);Ry=zeros(1,n1);
    for i=1:n1
        [~,r]=find(yd==ry(i));
        Ry(i)=r;% 获得y数据列秩统计量
    end
    xp=cumsum(ones(1,n1));
    d=xp-Ry;
    Qxy=1-6/(n1*(n1^2-1))*sum(d.^2);
    T=Qxy*sqrt(n1-2)/sqrt(1-Qxy^2);
    t_0=tinv(0.975,n1-2);% 计算上alpha/2分位数
    if(abs(T)<=t_0)
        break
    elseif(dp==dm+2)
        disp('原始数据列无法差分平稳化，可能不适合ARIMA模型')
        break
    else
        yd=diff(yd);
    end
end
d_p=dp-1;% 最终差分次数
% 自相关系数，偏自相关系数计算
y_mean=mean(yd);% 数据yd平均值
y_0=yd-y_mean;% 序列减去其均值
rho=zeros(1,n1);
for k=0:n1-1
    rho(k+1)=sum(y_0(1:n1-k).*y_0(k+1:n1))/sum(y_0.^2);
end
phi=zeros(1,n1-1);
for k=0:n1-2
    tic
    rho1=rho(2:k+2);
    kp=length(rho1);
    rho2=zeros(kp);
    for k0=1:kp
        rho2(k0,k0:end)=rho(1:kp-k0+1);
        if(k0>=2)
            rho2(k0,1:k0-1)=fliplr(rho1(1:k0-1));
        end
    end
    phi1=rho2\rho1';phi(k+1)=phi1(end);
    toc
end
% 定阶预测
% ARMA(p,q)预测一般需满足：数据自相关系数(ACF)q阶后衰减趋于0，偏自相关系数
% (PACF)p阶后衰减趋于0.
% 如果样本ACF和样本PACF在最初k阶明显大于2倍标准差，而后几乎95%的系数都落在2
% 倍标准差的范围内，且非零系数衰减为小值波动的过程非常突然，通常视为k阶截尾；
% 如果有超过5%的样本相关系数大于2倍标准差，或者非零系数衰减为小值波动的过程
% 比较缓慢或连续，通常视为拖尾。
rho_std=std(rho);% 自相关系数标准差
phi_std=std(phi);% 偏自相关系数标准差
p=0;q=0;
for i=rho
    if(abs(i)>2*rho_std)
        p=p+1;
    else
        break
    end
end
for i=phi
    if(abs(i)>2*phi_std)
        q=q+1;
    else
        break
    end
end
% 以上程序简单进行了定阶，通常利用贝叶斯信息等可以更准确的进行定阶，限于篇幅
% 本文不再赘述。
% 利用matlab内置函数进行预测；
AR_Order=p;
MA_Order=q;
AR=arima(AR_Order, d_p, MA_Order);
EstMdl=estimate(AR,y');
step=100;% 预测步数
fore_y=forecast(EstMdl,step,'Y0',y');
figure(1)
plot([y,fore_y'],'r*-')
hold on
plot(y,'b*-')
hold off
legend('预测值','实际值')
% 由结果分析可知，ARIMA模型适合短期的预测，且阶数的选择会影响预测结果。
%% 4.插值与拟合Interpolation and fitting
%%  4.1.插值
% 若数据有多项式特征，则可以用多项式函数进行插值。可利用待定系数法确定插值多
% 项式。
clear,clc
% 设有待插值原始数据列：
x0=(1:6)';y0=[16,18,21,17,15,12]';
A=vander(x0);% 返回x0范德蒙德矩阵
p=A\y0;% 求出多项式系数
x=0:0.1:6;
yh=polyval(p,x);plot(x,yh)% 绘制函数插值图像
% 用多项式做插值函数，随着插值节点的增加，插值多项式的次数越多，高次插值不但
% 计算复杂而且会产生龙格振荡现象。因此实际情况中一般会采用其他插值方法。
% 分段线性插值：将相邻节点用直线连接起来；
% 三次样条插值：插值函数在每个子区间均为三次多项式，且满足一定边界条件。
% 本节介绍matlab中一维插值函数：
clear,clc
% 设有待插值原始数据列：
x0=(1:6)';y0=[16,18,21,17,15,12]';
method='spline';% 设置插值方法，可以为'nearest'最近邻插值；'linear'线性插值；
% 'spline'三次样条插值；'cubic'立方插值
xq=0:0.1:6;% 待求点位置
vq=interp1(x0,y0,xq,method);plot(xq,vq)
% 二维插值函数interp2与interp1用法基本相同。
%%  4.2.拟合
% 已知一组二维数据(x(i),y(i)),i=1,2,...,n，要寻找一个函数y=f(x)使f(x)在某准
% 则下与所有数据点最为接近。此时可以采用偏差平方和最小准则。即：
% min J=sum((f(1:n)-y(1:n)).^2);这一原则称为最小二乘元则。
% 线性最小二乘：
% 给定一个线性无关的函数系{phi(k)|,k=1,2,...,m}，若拟合函数可以以其线性组合
% 形式出现:f(x)=sum(a(1:m).*phi(1:m))，则f(x)为关于系数{a(k)|,k=1,2,...,m}
% 的线性函数。
% 对于拟合函数：f(x)=a(1)*phi(1)+a(2)*phi(2)+,...,+a(m)*phi(m);
% 记：R=[phi(1:m)|x(1);phi(1:m)|x(2);...;phi(1:m)|x(n)];A=[a(1:m)];Y=[y(1:n)];
% 有：A=(R'*R)\R'*Y。
% matlab求约束线性最小二乘解：
clear,clc
x=[3,5,6,7,4,8,5,9]';y=[4,9,5,3,8,5,8,5]';% 观测数据
c=[exp(x),log(x)];% 线性无关函数系
a=lsqlin(c,y);% 拟合参数
% matlab求多项式拟合：
clear,clc
x=[3,5,6,7,4,8,5,9]';y=[4,9,5,3,8,5,8,5]';% 观测数据
p=polyfit(x,y,3);% 拟合3次多项式，返回值p为多项式系数，从高次幂到低次幂
% matlab fit和fittype函数
clear,clc
x=[3,5,6,7,4,8,5,9]';y=[4,9,5,3,8,5,8,5]';% 观测数据
g=fittype('a*exp(x)+b*log(x)','dependent',{'y'},'independent',{'x'});
f=fit(x,y,g,'StartPoint',rand(1,2));